home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48hor2 / fft7.src < prev    next >
Text File  |  1992-08-18  |  29KB  |  875 lines
  1. (Comp.sources.hp48)
  2. Item: 140 by kcooke@milton.u.washington.edu [Ken Cooke]
  3. Subj: Machine Code FFT source code
  4. Date: 14 Jul 1992
  5.  
  6. I'm glad to hear that someone liked the machine code FFT programs.
  7. In response to email requests, here is the SASM source code.  It can
  8. be assembled with the HP tools.  I didn't include FFT, because it is a
  9. subset of IFFT (I commented the lines to remove).
  10.  
  11. I apologize for the shortage of higher-level comments (most of this is
  12. still on paper).  Also, I traded some readability to avoid subroutine calls
  13. (hardware stack levels are very precious) and keep variables in registers.
  14. A few of the ROM calls are "unsupported entry points", which may need
  15. fixing in the future.  With rumors of new 48's in the works, ones with
  16. machine code Eqn Writers, re-written plotter user interfaces, and maybe
  17. even faster Saturn processors...well, lets just say that it may be wise to
  18. stick to entries listed in ENTRIES.A.
  19.  
  20. I hope this is of use to someone interested in ML math software.  While I
  21. like playing games on my 48 as much as the next guy, it would be nice to 
  22. see more math tools using the capabilities of sys-RPL and ML.  (I do, 
  23. occasionally, use my 48 as a calculating device :)  Unfortunately,
  24. there is no documentation on using the low-level floating point routines 
  25. (MULTF, RADDF, etc) so a little hacking with Voyager or MLDL is required.
  26.  
  27. I would love to hear from other math-hackers out there...
  28.  
  29. Ken Cooke
  30. N7VFE
  31. kcooke@u.washington.edu
  32.  
  33. ----------
  34.  
  35. * ANGLE.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
  36. * Copyright 1992 Ken Cooke
  37. * Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
  38. * same size, that contains the angle of each complex element,
  39. * according to the current trig mode.
  40.  
  41.          TITLE ANGLE
  42. ASSEMBLE
  43.          NIBASC    /HPHP48-A/
  44. RPL
  45. ::
  46. CK1NoBlame
  47. CK&DISPATCH1
  48. FOUR
  49. ::
  50. * IF REAL ARRAY, EXIT UNCHANGED
  51.          DUP TYPECARRY? NOT?SEMI                         ( ORIG )
  52. * MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
  53.          DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY         ( ORIG ANS )
  54. CODE
  55. ****************** UNLISTED ENTRY POINTS *****************
  56. ATTN?Lp  = #0CA88
  57. ARG15    = #2B6D7
  58. SetTrigMode        = #2AEF6
  59. ************************ MCODE ANGLE ***********************
  60. * ENTRY:
  61. *       STK2: ORIGINAL ARRY1 (COMPLEX)
  62. *        STK1: ANSWER ARRY (REAL)
  63.  
  64. * EXIT: STK1: ANSWER ARRAY (REAL)
  65.  
  66.          NIBHEX  823        ;CLEAR SB AND XM IN HST
  67.          D1=D1+    5        ;DROP STK1
  68.          D=D+1     A        ;RETURN MEM
  69.          GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK
  70.  
  71.          A=DAT1    A        ;A->ORIG
  72.          D1=D1-    5        ;POINT BACK TO ANS
  73.          C=DAT1    A        ;C->ANS
  74.          R3=C.F    A        ;SAVE ->ANS IN R3
  75.          D0=C
  76.          D0=D0+    15
  77.          D0=D0+    10       ;D0->ANS DATA (OR DIM2 IF ARRAY)
  78.  
  79.          D1=A               ;D1->ORIG
  80.          D1=D1+    15       ;D1->#DIMS
  81.          A=DAT1    A        ;A=#DIMS
  82.          D1=D1+    5        ;D1->DIM1
  83.          C=DAT1    A        ;C=DIM1
  84.          D1=D1+    5        ;D1->ORIG DATA (OR DIM2 IF ARRAY)
  85.          A=A-1     A        ;CHECK FOR VECTOR
  86.          ?A=0      A
  87.          GOYES     ISVECT
  88. ***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
  89.          A=DAT1    A        ;A=DIM2
  90.          D1=D1+    5        ;D1->ORIG DATA
  91.          D0=D0+    5        ;D0->ANS DATA
  92.          SETHEX
  93.          GOSBVL    =MUL#    ;B=A*C
  94.          C=B       A        ;C=DIM1*DIM2
  95.  
  96. ISVECT  C=C-1      A        ;SUB 1 SO CAN USE BORROW FOR LOOPING
  97.          R4=C.F    A        ;SAVE #ELEMENTS-1 IN R4
  98.  
  99. ABSLOOP GOSBVL     =ATTN?Lp ;ABORT IF ATTN (ON) WAS PRESSED
  100.          A=DAT1    W        ;A=RE
  101.          D1=D1+    16       ;D1->IM
  102.          CD0EX
  103.          RSTK=C             ;PUSH D0
  104.          C=DAT1    W        ;C=IM
  105.          D1=D1+    16       ;D1->NEXT RE
  106.          SETDEC
  107.          GOSBVL    =SPLTAC  ;CONVERT BOTH TO LONGS
  108.          GOSBVL    =SetTrigMode       ;SET ST ACCORDING TO TRIG MODE FLAGS
  109.          GOSBVL    =ARG15   ;AB=ANGLE (TRASHES A,B,C,D,D0,R0,R1)
  110.          GOSBVL    =PACK    ;A=ABSVAL
  111.          C=RSTK
  112.          D0=C               ;POP D0
  113.          DAT0=A    W        ;STORE ABSVAL IN ANSWER ARRAY
  114.          D0=D0+    16       ;POINT AT NEXT ANSWER
  115. * CHECK FOR DONE
  116.          C=R4.F    A        ;GET COUNTER
  117.          SETHEX
  118.          C=C-1     A
  119.          R4=C.F    A        ;DEC COUNTER
  120.          GONC    ABSLOOP    ;LOOP UNTIL BORROW
  121.  
  122. ***NOW ANSWER CONTAINS ABS VALUES OF ORIG
  123.          A=R3.F    A
  124.          GOSBVL    =GPOverWrALp       ;GETPTR, REPLACE TOS WITH ANS, LOOP
  125.  
  126. ENDCODE
  127. ;
  128. ;
  129. END_SRC
  130.  
  131.  
  132. BEGIN_SRC ifft.s
  133. ******************************************************************
  134. * MACHINE CODE IFFT by Ken Cooke
  135. * Copyright 1992
  136. ******************************************************************
  137.  
  138. * RADIX-2 DIT IFFT
  139. * Input: real or complex vector on STK1, with size=2^n
  140. * Ouput: complex vector containing IFFT
  141. * Note: uses the definition that FFT uses exp(-xxx) and IFFT uses exp(+xxx)
  142.  
  143.          TITLE IFFT
  144. ASSEMBLE
  145.          NIBASC    /HPHP48-A/
  146. RPL
  147. ::
  148. CK1NoBlame
  149. CK&DISPATCH1
  150. FOUR
  151. ::
  152.          DUP MDIMS                              ( ARRY #N flag )
  153. * CHECK FOR MATRIX                              (NcaseDIMERR = 37DF6 )
  154.          NOT NcaseDIMERR            ( ARRY #N )
  155. * CHECK FOR SIZE=2^N
  156.          DUP                                    ( ARRY #N #N )
  157.  
  158. CODE
  159. ****************** CHECK THAT SIZE=POWER OF TWO **************
  160. *CHKSIZE2 REPLACES BINT (ON TOS) WITH TRUE IF BINT IS A
  161. * POWER OF 2, FALSE OTHERWISE.
  162. *NOTE: WORKS UP TO 2^19, AND RETURNS TRUE FOR BINT=0
  163. WrTLoop  EQU       #03B1A
  164. WrFLoop  EQU       #03B06
  165.          C=DAT1    A        ;C->BINT
  166.          CD1EX              ;SAVE D1 IN C.A
  167.          D1=D1+    5        ;D1->DATA
  168.          A=DAT1    A        ;A.A=BINT
  169.          D1=C               ;RESTORE D1
  170.  
  171.          SETHEX
  172.          C=A       A
  173.          A=-A-1    A        ;A=ONES COMP
  174.          C=-C      A        ;C=TWOS COMP
  175.          A=A!C     A        ;NO ZEROS IF BINT=2^N
  176.          A=A+1     A        ;WILL SET CARRY IF BINT=2^N
  177.          GOC       DoTrue
  178.          GOVLNG    WrFLoop  ;OVERWRITE BINT WITH FALSE AND CONT
  179. DoTrue   GOVLNG    WrTLoop  ;OVERWRITE BINT WITH TRUE AND CONT
  180. ENDCODE
  181.  
  182. *                  ( ARRY #N FLAG )
  183. NcaseDIMERR        ( IF FALSE THEN "INVALID DIMENSION" ERROR )
  184.  
  185. * IF SIZE=1, RETURN WITH ARRY ON STACK
  186.          DUP#1= ITE DROP    ( ARRY #N )
  187.  
  188. ::
  189. * CREATE TEMP STORAGE (MUST BE 172 NIBS BIGGER THAN UNPACKED A FOR VARS)
  190.  
  191.          FORTYTWO OVER #* 172 #+ MAKEHEX$ SWAP  ( ARRY HEX$ #N )
  192.          ONE{}N C%0 MAKEARRY                             ( ARRY HEX$ ANSARRY )
  193.  
  194. * DO THE FFT
  195.  
  196. CODE
  197. ****************** UNLISTED ENTRY POINTS *****************
  198. PUTAB1   = #2C066
  199. GETCD1   = #2C017
  200. aSIN15   = #2B6E0
  201. ATTN?Lp  = #0CA88
  202. longpi   = #2A45D
  203. varnibs  = 172
  204. ************************ MCODE FFT ***********************
  205. * ENTRY:
  206. *        STK3: (orig) real/complex data array
  207. *       STK2: (temp) hex$ (to hold vars and unpacked reals during FFT)
  208. *        STK1: (ans) complex answer array
  209.  
  210. * EXIT: STK1: complex FFT array (ans)
  211.  
  212.          NIBHEX  823        ;CLEAR SB AND XM IN HST
  213.          D1=D1+    10       ;DROP temp AND ans
  214.          D=D+1     A        ;RETURN MEM
  215.          D=D+1     A
  216.          GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK
  217.  
  218. * SAVE NEEDED POINTERS IN SCRATCH REGS
  219.          D1=D1-    5        ;POINT BACK TO temp
  220.          A=DAT1    A        ;A->temp
  221.          A=A+CON   A,10     ;SKIP HEADER TO vars
  222.          R4=A.F    A        ;R4->vars
  223.          P=        0        ;*****REMEMBER:ALWAYS HAVE P=0 BEFORE LC(5)*****
  224.          LC(5)     varnibs  ;
  225.          SETHEX
  226.          A=A+C     A        ;SKIP OVER VARS
  227.          R1=A.F    A        ;R1->data
  228.  
  229.          D1=D1+    5        ;POINT TO orig
  230.          A=DAT1    A        ;A->orig
  231.          D0=A
  232.          D0=D0+    10       ;D0->OBJ PROLOG
  233.  
  234. * SET FLAG IF REAL, FOR USE WHEN UNPACKING
  235. * USE D.15 INSTEAD OF ST TO AVOID TROUBLE
  236.          D=0       S        ;USE D.15 AS FLAG (0=COMPLEX,1=REAL)
  237.          A=DAT0    A        ;A=OB PROLOG
  238.          LC(5)     #02977   ;C=COMPLEX PROLOG
  239.          ?A=C      A
  240.          GOYES     CMPLEX
  241.          D=D+1     S        ;IF REAL, D.15=1
  242.  
  243. CMPLEX   D0=D0+    10       ;D0->SIZE
  244.          A=DAT0    A        ;A.A=SIZE
  245.          R2=A.F    A        ;STORE SIZE IN R2
  246.          D0=D0+    5        ;D0->DATA
  247.          CD0EX
  248.          R0=C.F    A        ;R0->origdata
  249.  
  250. *FIND LOG2 OF SIZE
  251.          C=0       S        ;C.S WILL HOLD SIZE2=LOGbase2(SIZE)
  252. LOG2LP   ASRB.F    A
  253.          C=C+1     S
  254.          ?ABIT=0   0
  255.          GOYES     LOG2LP
  256.          R2=C.F    S        ;STORE SIZE2 IN R2.15
  257.  
  258. *************************** BITREV AND UNPACK ****************************
  259. *NOW, CONVERT REALS IN orig INTO EXT.REALS IN data
  260. * AT THE SAME TIME, REARRANGE INTO BIT-REVERSED ORDER
  261.  
  262.          A=R0.F    A
  263.          D0=A               ;D0->origdata
  264.          D=0       A        ;D HOLDS INDEX INTO orig (START AT 0)
  265.  
  266. BITREV   CPEX      15       ;PUT BIT-COUNTER IN P
  267.          P=P-1           ;SUB ONE SO CAN USE BORROW FOR BITLOOP
  268.          C=D       A        ;C GETS INDEX INTO orig
  269.          B=0       A        ;B WILL BE BIT-REV INDEX
  270.  
  271. *BIT REVERSAL
  272. BITLOOP B=B+B      A        ;SHIFT B LEFT ONE BIT
  273.          ?CBIT=0   0        ;TEST LSB
  274.          GOYES     NOCARRY  ;SKIP IF ZERO
  275.          B=B+1     A        ;ELSE FEED INTO B
  276. NOCARRY  CSRB.F    A       ;SHIFT C RIGHT ONE BIT
  277.          P=P-1              ;DEC BIT-COUNTER (BORROW WHEN DONE)
  278.          GONC      BITLOOP  ;LOOP IF MORE
  279.  
  280. *MULT B (BITREV INDEX) BY 21 AND ADD TO ->data
  281.          A=B       A        ;COULD USE =MUL# BUT THIS IS FASTER/SHORTER
  282.          B=B+B     A
  283.          B=B+B     A        ;B=4B
  284.          A=A+B     A        ;A=5B
  285.          B=B+B     A
  286.          B=B+B     A        ;B=16B
  287.          B=A+B     A        ;B=21B
  288.          B=B+B     A        ;B=42B
  289.  
  290.          A=R1.F    A        ;A->data
  291.          A=A+B     A        ;
  292.          D1=A               ;D1 POINTS INTO data AT BITREV POSITION
  293.          A=DAT0    W        ;GET RE(Ai)
  294.          D0=D0+    16       ;POINT AT NEXT REAL
  295.          GOSBVL    =SPLITA  ;NO NEED TO SETDEC (DONE IN SPLITA)
  296.          GOSBVL  =PUTAB1    ;PUT LONG AT DAT1, D1+=21
  297.  
  298. * IF INPUT ARRAY WAS REAL, USE IMAG=0
  299.          ?D=0      S        ;IF COMPLEX
  300.          GOYES     GETIM    ; GET IMAG PART
  301.          A=0       W        ;ELSE USE IMAG=0
  302.          B=0       W
  303.          GONC      NOGETIM  ;GO ALWAYS
  304.  
  305. GETIM    A=DAT0    W        ;GET IM(Ai)
  306.          D0=D0+    16
  307.          GOSBVL    =SPLITA
  308.  
  309. NOGETIM  GOSBVL  =PUTAB1    ;STORE IMAG PART
  310.          SETHEX
  311.          D=D+1     A        ;INC INDEX
  312.          C=R2               ;LOAD SIZE AND SIZE2 INTO C
  313.          ?D=C      A        ;COMPARE INDEX TO SIZE
  314.          GOYES     BITDONE
  315.          GOTO      BITREV   ;REPEAT IF NOT DONE
  316. BITDONE
  317.  
  318. ********** data NOW HOLDS LONG BITREV ARRAY **********
  319.  
  320.  
  321. ******************** BUTTERFLIES *******************
  322. * ENTRY: temp CONTAINS BIT-REV LONG COMPLEX DATA
  323. *        A IS NO LONGER NEEDED
  324. * Define template for var storage: (var names from Num.Recipes)
  325.  
  326. istep    = 0
  327. n        = (istep)+5
  328. mmax     = (n)+5
  329. data0    = (mmax)+5
  330. wr       = (data0)+5
  331. wi       = (wr)+21
  332. tempr    = (wi)+21
  333. wtemp    = (tempr)+21
  334. theta    = (wtemp)+21
  335. wpi      = (theta)+21
  336. wpr      = (wpi)+21
  337. m        = (wpr)+21
  338.  
  339. *FIRST, STORE SCRATCH REGS IN VARS
  340.  
  341.          SETHEX
  342.          A=R4.F    A        ;A->istep
  343.          D1=A
  344.          D1=D1+    5        ;D1->n
  345.          C=R2.F    A        ;C=size
  346.          C=C+C     A        ;C=2*size
  347.          DAT1=C    A        ;n=2*size
  348.  
  349.          D1=D1+    5        ;D1->mmax
  350.          P=        0
  351.          LC(5)     #2       ;C=2
  352.          DAT1=C    A        ;mmax=2
  353.  
  354.          D1=D1+    5        ;D1->data0          ***STORE POINTER TO data[0] HERE
  355.          C=R1.F    A        ;***THIS CORRECTS FOR "FORTRAN INDEXING"
  356.          C=C-CON   A,16     ;*** SINCE ALGORITHM ASSUMES INDEXES START AT 1
  357.          C=C-CON   A,5      ;*** I.E. data[1] = data0+21 = FIRST SPOT
  358.          DAT1=C    A
  359.  
  360. **REDUCED TO ONE SINE/MAINLOOP BY USING wtemp TO HOLD LAST SIN(theta)
  361. ** AND CALCULATING SIN(theta/2)
  362.          LC(5)     wtemp   ;INIT wtemp TO AVOID SIN(PI) SINCE NOT RETURN ZERO
  363.          C=C+A     A        ;A STILL =R4
  364.          D0=C               ;D0->wtemp
  365.          A=0       W
  366.          B=0       W        ;AB=0.0
  367.          GOSBVL    =PUTAB0  ;STORE wtemp=0.0, D0->theta
  368.  
  369. **INIT theta TO -pi FOR FFT, +pi FOR IFFT
  370.          D1=(5)  =longpi
  371.          GOSBVL    =GETAB1  ;get pi from ROM
  372. *        SETDEC                                 *********enable for FFT
  373. *        A=-A-1    S        ;A=-PI              *********enable for FFT
  374.          GOSBVL    =PUTAB0  ;STORE IN theta
  375.  
  376. ************************ BEGIN MAIN LOOP ***********************
  377. MAINLOOP
  378.          C=R4.F    A        ;C->vars
  379.          D1=C
  380.          D1=D1+    5        ;D1->n
  381.          C=DAT1    A        ;C=n
  382.          D1=D1+    5        ;POINT AT mmax
  383.          A=DAT1    A        ;A=mmax
  384. * MAINLOOP TEST
  385.          ?A<C      A        ;IF (mmax>=n)
  386.          GOYES     MAINCONT ; THEN DONE
  387.          GOTO      MAINDONE
  388. ** A=mmax, D1->mmax
  389. MAINCONT
  390.          SETHEX
  391.          A=A+A     A        ;A=2*mmax
  392.          D1=D1-    10       ;istep
  393.          DAT1=A    A        ;STORE istep=2*mmax
  394.  
  395.          D1=D1+    10
  396.          D1=D1+    10       ;D1->wr
  397.          A=0       W
  398.          B=0       W
  399.          P=        14       ;CREATE AB=1.0
  400.          B=B+1     P
  401.          GOSBVL    =PUTAB1  ;STORE wr=1.0, D1->wi
  402.          B=0       P
  403.          GOSBVL    =PUTAB1  ;STORE wi=0.0, D1->tempr
  404.  
  405.          D1=D1+    16
  406.          D1=D1+    5        ;D1->wtemp
  407.          CD1EX
  408.          RSTK=C             ;PUSH ->wtemp
  409.          D0=C               ;D0->wtemp
  410.          D1=C               ;D1->wtemp
  411.          D1=D1+    16
  412.          D1=D1+    16
  413.          D1=D1+    10       ;D1->wpi
  414.  
  415.          SETDEC             :***NOTE: SETDEC BEFORE ALL FP MATH (WILL STAY IN DEC)
  416.          GOSBVL    =GETAB0  ;AB=wtemp, D0->theta
  417.          GOSBVL    =PUTAB1  ;STORE wpi=wtemp, D1->wpr
  418.          GOSBVL    =GETAB0  ;AB=theta, D0->wpi
  419.          GOSBVL    =DIV2    ;AB=theta/2
  420.          D0=D0-    16
  421.          D0=D0-    5        ;D0->theta
  422.          GOSBVL    =PUTAB0  ;STORE theta=theta/2
  423.          ST=1      9        ;SET UP FOR SIN (IN RADIANS)
  424.          ST=0      4
  425.          GOSBVL    =aSIN15  ;**NOTE:TRASHES A,B,C,D,P,D0,R0,R1!!!
  426.          C=RSTK             ;POP ->wtemp
  427.          D0=C               ;D0->wtemp
  428.          GOSBVL    =PUTAB0  ;STORE wtemp=SIN(theta)
  429.  
  430.          C=B       W
  431.          D=C       W        ;CD=AB
  432.          C=A       W
  433.          GOSBVL    =MULTF   ;SQUARE wtemp
  434.          C=B       W
  435.          D=C       W        ;CD=AB
  436.          C=A       W
  437.          GOSBVL    =RADDF   ;DOUBLE AB (D0 TRASHED)
  438.          A=-A-1    S        ;NEGATE AB
  439.          GOSBVL    =PUTAB1  ;STORE wpr=-2*SQR(wtemp), D1->m
  440.  
  441.          A=0       A
  442.          A=A+1     B        ;A=1
  443.          DAT1=A    A        ;STORE m=1
  444.  
  445. ************************ BEGIN OUTER LOOP ***********************
  446. OUTLOOP  SETHEX
  447.          A=R4.F    A
  448.          D0=A
  449.          D0=D0+    10       ;D0->mmax
  450.          P=        0
  451.          LC(5)     m
  452.          C=C+A     A
  453.          D1=C               ;D1->m
  454.          A=DAT0    A        ;A=mmax
  455.          C=DAT1    A        ;C=m
  456. *OUTLOOP TEST
  457.          ?C<A      A        ;IF (m>=mmax)
  458.          GOYES     NOTDONE  ; THEN DONE
  459.          GOTO      OUTDONE
  460. NOTDONE  RSTK=C             ;PUSH i=m (i KEPT ON TOS)
  461.          A=R4.F    A
  462.          D1=A               ;D1->n-5
  463.  
  464. ************************ BEGIN INNER LOOP ***********************
  465. INLOOP
  466.          GOSBVL    =ATTN?Lp ;IF ATTN(ON) HAS BEEN PRESSED, GETPTR AND LOOP
  467.          D1=D1+    5        ;D1->n
  468.          A=DAT1    A        ;A=n
  469.          C=RSTK             ;C=i
  470. *INLOOP TEST
  471.          ?C<=A     A        ;IF (i>n)
  472.          GOYES     INCONT   ; THEN DONE
  473.          GOTO      INDONE
  474. INCONT   RSTK=C             ;PUSH i BACK
  475.          D1=D1+    5        ;D1->mmax
  476.          A=DAT1    A        ;A=mmax
  477.          SETHEX
  478.          A=A+C     A        ;A=j
  479.          B=A       A
  480.          B=B+B     A
  481.          B=B+B     A
  482.          A=A+B     A
  483.          B=B+B     A
  484.          B=B+B     A
  485.          A=A+B     A        ;A=21*j
  486.  
  487.          B=C       A
  488.          C=C+C     A
  489.          C=C+C     A
  490.          B=B+C     A
  491.          C=C+C     A
  492.          C=C+C     A
  493.          B=C+B     A        ;B=21*i
  494.  
  495.          D1=D1+    5        ;D1->data0
  496.          C=DAT1    A        ;C->data
  497.          A=A+C     A        ;A->data[j]
  498.          C=B+C     A        ;C->data[i]
  499.          RSTK=C             ;PUSH ->data[i]
  500.          C=A       A
  501.          RSTK=C             ;PUSH ->data[j]
  502.  
  503.          D1=D1+    5        ;D1->wr
  504.          D0=C               ;D0->data[j]
  505.          GOSBVL    =GETAB1  ;AB=wr, D1->wi
  506.          GOSBVL    =STAB2   ;SAVE wr IN R2/R3
  507.          GOSBVL    =GETCD0  ;CD=data[j], D0->data[j+1]
  508.          SETDEC
  509.          GOSBVL    =MULTF   ;AB=wr*data[j]
  510.          GOSBVL    =STAB0   ;SAVE IN R0/R1
  511.  
  512.          GOSBVL    =GETAB1  ;AB=wi, D1->tempr
  513.          GOSBVL    =GETCD0  ;CD=data[j+1]
  514.          GOSBVL    =MULTF   ;AB=wi*data[j+1]
  515.          GOSBVL    =RCCD0   ;CD=wr*data[j]
  516.          A=-A-1    S        ;NEGATE AB
  517.          GOSBVL    =RADDF   ;AB=RESULT (D0 TRASHED)
  518.  
  519.          GOSBVL    =PUTAB1  ;tempr=RESULT, D1->tempi
  520.          C=RSTK
  521.          RSTK=C             ;COPY ->data[j] FROM TOS
  522.          D0=C               ;D0->data[j]
  523.          GOSBVL    =GETCD0  ;CD=data[j], D0->data[j+1]
  524.          D1=D1-    16
  525.          D1=D1-    16
  526.          D1=D1-    10       ;D1->wi
  527.          GOSBVL    =GETAB1  ;AB=wi
  528.  
  529.          GOSBVL    =MULTF
  530.          GOSBVL    =STAB0   ;STORE wi*data[j] IN R0/R1
  531.          GOSBVL    =GETCD0  ;CD=data[j+1]
  532.          GOSBVL    =RCAB2   ;AB=wr
  533.          GOSBVL    =MULTF   ;AB=wr*data[j+1]
  534.          GOSBVL    =RCCD0   ;CD=wi*data[j]
  535.          GOSBVL    =RADDF   ;AB=tempi (D0 TRASHED)
  536.  
  537.          GOSBVL    =STAB0   ;SAVE tempi IN R0/R1
  538.          C=RSTK
  539. ***STORE D1=>data[j]
  540.          D1=C               ;D1->data[j]
  541.          C=RSTK
  542.          RSTK=C             ;COPY ->data[i] FROM TOS
  543.          C=C+CON   A,16
  544.          C=C+CON   A,5      ;C->data[i+1]
  545.          RSTK=C             ;PUSH ->data[i+1]
  546.          D0=C
  547.          GOSBVL    =GETCD0  ;CD=data[i+1]
  548.          GOSBVL    =STCD2   ;SAVE data[i+1] IN R2/R3
  549.  
  550.          A=-A-1    S        ;NEGATE tempi
  551.          GOSBVL    =RADDF   ;AB=data[i+1] - tempi
  552.          CD1EX
  553.          D1=C               ;C->data[j]
  554.          D0=C
  555.          D0=D0+    16
  556.          D0=D0+    5        ;D0->data[j+1]
  557.  
  558.          GOSBVL    =PUTAB0  ;STORE data[j+1]
  559.          GOSBVL    =RCAB0   ;AB=tempi
  560.          GOSBVL    =RCCD2   ;CD=data[i+1]
  561.          GOSBVL    =RADDF   ;AB=tempi+(data[i+1])
  562.          C=RSTK             ;POP ->data[i+1]
  563.          D0=C               ;D0->data[i+1]
  564.          GOSBVL    =PUTAB0  ;STORE data[i+1]
  565.  
  566.          A=R4.F    A
  567.          P=        0
  568.          LC(5)     tempr
  569.          SETHEX
  570.          A=A+C     A
  571.          D0=A               ;D0->tempr
  572.          GOSBVL    =GETAB0  ;AB=tempr
  573.          GOSBVL    =STAB0   ;STORE tempr IN R0/R1
  574.          C=RSTK
  575.          RSTK=C             ;COPY ->data[i] FROM TOS
  576.          D0=C               ;D0->data[i]
  577.          GOSBVL    =GETCD0  ;CD=data[i]
  578.          GOSBVL    =STCD2   ;SAVE data[i] IN R2/R3
  579.          SETDEC
  580.          A=-A-1    S        ;NEGATE tempr
  581.          GOSBVL    =RADDF   ;AB=data[i]-tempr
  582. **USE D1=>data[j]
  583.          GOSBVL    =PUTAB1  ;STORE data[j]
  584.  
  585.          GOSBVL    =RCAB0   ;AB=tempr
  586.          GOSBVL    =RCCD2   ;CD=data[i]
  587.          GOSBVL    =RADDF   ;AB=data[i]+tempr
  588.          C=RSTK             ;POP ->data[i]
  589.          D0=C               ;D0->data[i]
  590.          GOSBVL    =PUTAB0  ;STORE data[i]
  591.  
  592.          A=R4.F    A
  593.          D1=A               ;D1->istep
  594.          A=DAT1    A        ;A=istep
  595.          C=RSTK             ;POP i
  596.          SETHEX
  597.          C=C+A     A        ;i=i+istep
  598.          RSTK=C             ;PUSH NEW i
  599.  
  600.          GOTO      INLOOP   ;NEXT ITERATION
  601. ************************ END INNER LOOP ***********************
  602. INDONE
  603.          A=R4.F    A
  604.          P=        0
  605.          LC(5)     wr
  606.          SETHEX
  607.          C=C+A     A
  608.          RSTK=C             ;PUSH *wr
  609.          D0=C               ;D0->wr
  610.          LC(5)     wpr
  611.          C=C+A     A
  612.          D1=C               ;D1->wpr
  613.  
  614.          SETDEC
  615.          GOSBVL    =GETAB0  ;AB=wr
  616.          GOSBVL    =STAB0   ;SAVE OLDwr IN R0/R1
  617.          GOSBVL    =GETCD1  ;CD=wpr, D1->wpr+21
  618.          GOSBVL    =MULTF   ;AB=wr*wpr
  619.          GOSBVL    =STAB2   ;SAVE IN R2/R3
  620.          GOSBVL    =GETCD0  ;CD=wi
  621.          D1=D1-    16
  622.          D1=D1-    16
  623.          D1=D1-    10       ;D1->wpi
  624.          GOSBVL    =GETAB1  ;AB=wpi, D1->wpr
  625.          GOSBVL    =MULTF   ;AB=wi*wpi
  626.          A=-A-1    S        ;NEGATE AB
  627.          GOSBVL    =RCCD2   ;CD=wr*wpr
  628.          GOSBVL    =RADDF   ;AB=wr*wpr - wi*wpi (D0 TRASHED)
  629.          GOSBVL    =RCCD0   ;CD=OLDwr
  630.          GOSBVL    =RADDF   ;AB=RESULT (D0 TRASHED)
  631.  
  632.          C=RSTK             ;POP *wr
  633.          D0=C               ;D0->wr
  634.          GOSBVL    =PUTAB0  ;STORE NEW wr, D0->wi
  635.          D1=D1-    16
  636.          D1=D1-    5        ;D1->wpi
  637.          GOSBVL    =GETAB1  ;AB=wpi, D1->wpr
  638.          GOSBVL    =RCCD0   ;CD=OLDwr
  639.          GOSBVL    =MULTF   ;AB=OLDwr*wpi
  640.          GOSBVL    =STAB0   ;SAVE IN R0/R1
  641.  
  642.          GOSBVL    =GETCD0  ;CD=wi
  643.          GOSBVL    =STCD2   ;SAVE wi IN R2/R3
  644.          GOSBVL    =GETAB1  ;AB=wpr, D1->m
  645.          GOSBVL    =MULTF   ;AB=wpr*wi
  646.          GOSBVL    =RCCD0   ;CD=OLDwr*wpi
  647.          GOSBVL    =RADDF   ;AB=wpr*wi + OLDwr*wpi (D0 TRASHED)
  648.          GOSBVL    =RCCD2   ;CD=wi
  649.          GOSBVL    =RADDF   ;AB=RESULT
  650.  
  651.          C=R4.F    A
  652.          D0=C
  653.          D0=D0+    16
  654.          D0=D0+    16
  655.          D0=D0+    9        ;D0->wi
  656.          GOSBVL    =PUTAB0  ;STORE NEW wi
  657.          C=DAT1    A        ;C=m
  658.          SETHEX
  659.          C=C+1     A
  660.          C=C+1     A        ;m=m+2
  661.          DAT1=C    A        ;STORE NEW m
  662.  
  663.          GOTO      OUTLOOP  ;NEXT ITERATION
  664. ************************ END OUTER LOOP ***********************
  665.  
  666. OUTDONE  SETHEX             ;D0 STILL POINTS TO mmax, A=mmax
  667.          A=A+A     A        ;A=2*mmax
  668.          DAT0=A    A        ;STORE NEW mmax
  669.  
  670.          GOTO      MAINLOOP
  671. ************************ END MAIN LOOP ***********************
  672.  
  673. MAINDONE
  674. * temp NOW HOLDS vars (FIRST 172 NIBS) AND LONG COMPLEX ARRAY OF ANSWERS
  675.  
  676. ***************** PACK AND COPY INTO ans *****************
  677.  
  678.          A=R4.F    A
  679.          D1=A               ;D1->vars
  680.          D1=D1+    5        ;D1->n
  681.          A=DAT1    A        ;A=n
  682.          R2=A.F    A        ;SAVE n IN R2 (COUNTER IN PACKLOOP)
  683.  
  684. ***REMOVE FOR FFT
  685.          ASRB.F    A        ;A=size
  686. ************* CONVERT A.A INTO EXT.REAL ************
  687. * USING HORNER'S RULE, 10100 = 2^4 + 2^2 = ((((1)2+0)2+1)2+0)2+0
  688. *  SINCE THE *2 AND ADDITION IS DONE IN DEC, RESULT IN DEC
  689. * ENTER: BINT IN A.A (ANY)
  690. * EXIT: EXT.REAL IN A/B (DEC)
  691. * TRASHES: A,B,C,P
  692.  
  693.          P=        0
  694.          LCHEX  19          ;REPEAT FOR EACH BIT IN A.A
  695.          B=0    W  ;B WILL HOLD MANTISSA
  696.  
  697. BITLP    B=B+B  W  ;B*=2
  698.          SETHEX
  699.          A=A+A  A  ;SHIFT MSB OF BINT INTO CARRY
  700.          SETDEC
  701.          GONC   SKIP
  702.          B=B+1  W
  703. SKIP     C=C-1  B  ;DEC BIT COUNTER
  704.          GONC   BITLP
  705.  
  706.          P=     5  ;A.A IS 0 FROM SHIFTING IN ZEROS
  707.          A=A-1  X
  708. NORMLP   A=A+1  X  ;A ACCUM'S EXPONENT BY COUNTING DIGITS
  709.          BSRC      ;ROTATE UNTIL MANT IS LEFT-JUSTIFIED
  710.          ?B#0   WP
  711.          GOYES  NORMLP
  712.  
  713.          BSR    W  ;SHIFT MANT INTO PLACE (B.S=0)
  714.          A=0    S  ;SIGN IS POSITIVE
  715. * AB NOW HOLDS EXT.REAL
  716. **************** END OF CONVERSION ****************
  717.          GOSBVL    =STAB0   ;SAVE SIZE IN R0/R1
  718. *******REMOVE FOR FFT
  719.  
  720.          SETHEX
  721.          A=R4.F    A        ;**THIS IS THE LAST TIME WE USE R4(temp)
  722.         P=         0
  723.          LC(5)     varnibs  ;SKIP OVER VARS TO DATA
  724.          C=A+C     A        ;C->data
  725.          D0=C               ;D0->data (SOURCE OF COPY)
  726.  
  727.          D1=(5)    =DSKTOP  ;RECALL SAVED D1
  728.          C=DAT1    A
  729.          D1=C               ;D1->STK3 (orig)
  730.          D1=D1-    10       ;D1->STK1 (ans)
  731.          C=DAT1    A
  732.          R3=C.F    A        ;SAVE ->ans IN R3
  733.          D1=C               ;D1->ans
  734.          D1=D1+    15      ;SKIP HEADER
  735.          D1=D1+    10       ;D1->ansdata (DEST OF COPY)
  736.  
  737. * NOW PACK AND COPY
  738. * NOTE: OVER/UNDERFLOW ERRORS ARE TRAPPED IN =PACK. orig IS SAVED AS TOS
  739. PAKLOOP
  740.          GOSBVL    =GETAB0  ;GET NEXT LONG
  741.          GOSBVL    =RCCD0   ;RECALL SIZE                 *REMOVE FOR FFT
  742.          SETDEC
  743.          GOSBVL    =DIVF    ;DIVIDE BY SIZE              *REMOVE FOR FFT
  744.          GOSBVL    =PACK    ;EXT.REAL->REAL
  745.  
  746.          DAT1=A    W        ;PUT REAL IN PLACE
  747.          D1=D1+    16       ;POINT AT NEXT REAL
  748.  
  749.          C=R2.F    A        ;GET COUNTER
  750.          SETHEX
  751.          C=C-1     A        ;DEC COUNTER
  752.          R2=C.F    A
  753.          ?C#0      A
  754.          GOYES     PAKLOOP
  755.  
  756.  
  757. * NOW RESULT IS PACKED INTO ans ARRAY
  758. * AT LAST, REPLACE STK1 (ORIGINAL ARRAY) WITH RESULT ARRAY
  759.          A=R3.F    A                  ;A->ans
  760.          GOVLNG    =GPOverWrALp       ;RESTORE REGS, PUT @A ON TOS, CONT RPL
  761. ENDCODE
  762. ;
  763. ;
  764. ;
  765. END_SRC
  766.  
  767.  
  768.  
  769. BEGIN_SRC absv.s
  770. * ABSV.S  WORKS FOR ANY SIZE VECTOR OR ARRAY
  771. * Copyright 1992 Ken Cooke
  772. * Takes a VECTOR/ARRAY in STK1.  Returns a real VECTOR/ARRAY of
  773. * same size, that contains the element-by-element absolute values.
  774.  
  775.          TITLE ABSV
  776. ASSEMBLE
  777.          NIBASC    /HPHP48-A/
  778. RPL
  779. ::
  780. CK1NoBlame
  781. CK&DISPATCH1
  782. FOUR
  783. ::
  784. * IF REAL ARRAY, EXIT UNCHANGED
  785.          DUP TYPECARRY? NOT?SEMI                         ( ORIG )
  786. * MALLOC ANSWER ARRAY (SAME SIZE AS ORIG BUT REAL )
  787.          DUP MDIMS ITE TWO{}N ONE{}N %0 MAKEARRY         ( ORIG ANS )
  788. CODE
  789. ****************** UNLISTED ENTRY POINTS *****************
  790. ATTN?Lp  = #0CA88
  791. SQR15    = #2B9F3
  792.  
  793. ************************ MCODE ABSV ***********************
  794. * ENTRY:
  795. *       STK2: ORIGINAL ARRY1 (COMPLEX)
  796. *        STK1: ANSWER ARRY (REAL)
  797.  
  798. * EXIT: STK1: ANSWER ARRAY (REAL)
  799.  
  800.          NIBHEX  823        ;CLEAR SB AND XM IN HST
  801.          D1=D1+    5        ;DROP STK1
  802.          D=D+1     A        ;RETURN MEM
  803.          GOSBVL    =SAVPTR  ;SAVE REGS WITH ONLY orig ON STACK
  804.  
  805. * NOW IF USER HITS "ON" KEY, WILL LEAVE ORIGINAL STACK
  806.          A=DAT1    A        ;A->ORIG
  807.          D1=D1-    5        ;POINT BACK TO ANS
  808.          C=DAT1    A        ;C->ANS
  809.          R0=C.F    A        ;SAVE ->ANS IN R0
  810.          D0=C
  811.          D0=D0+    15
  812.          D0=D0+    10       ;D0->ANS DATA (OR DIM2 IF ARRAY)
  813.  
  814.          D1=A               ;D1->ORIG
  815.          D1=D1+    15       ;D1->#DIMS
  816.          A=DAT1    A        ;A=#DIMS
  817.          D1=D1+    5        ;D1->DIM1
  818.          C=DAT1    A        ;C=DIM1
  819.          D1=D1+    5        ;D1->ORIG DATA (OR DIM2 IF ARRAY)
  820.          A=A-1     A        ;CHECK FOR VECTOR
  821.          ?A=0      A
  822.          GOYES     ISVECT
  823. ***IF ARRAY, #ELEMENTS=DIM1*DIM2.  ALSO, SKIP OVER DIM2 FIELD
  824.          A=DAT1    A        ;A=DIM2
  825.          D1=D1+    5        ;D1->ORIG DATA
  826.          D0=D0+    5        ;D0->ANS DATA
  827.          SETHEX
  828.          GOSBVL    =MUL#    ;B=A*C
  829.          C=B       A        ;C=DIM1*DIM2
  830.  
  831. ISVECT  C=C-1      A        ;SUB 1 SO CAN USE BORROW FOR LOOPING
  832.          R4=C.F    A        ;SAVE #ELEMENTS-1 IN R4
  833.  
  834. ABSLOOP GOSBVL     =ATTN?Lp ;ABORT IF ATTN (ON) WAS PRESSED
  835.          SETDEC
  836.          GOSUB     SQDAT1   ;SQUARE RE PART
  837.          GOSBVL    =STAB2   ;SAVE IN R2/R3
  838.          GOSUB     SQDAT1   ;SQUARE IM PART
  839.  
  840.          CD0EX
  841.          RSTK=C             ;PUSH D0
  842.          GOSBVL    =RCCD2   ;CD=SQR(RE)
  843.          GOSBVL    =RADDF   ;AB=SQR(RE)+SQR(IM) (D0 TRASHED)
  844.          GOSBVL    =SQR15   ;AB=LONG ABSVAL
  845.          GOSBVL    =PACK    ;A=ABSVAL
  846.          C=RSTK
  847.          D0=C               ;POP D0
  848.          DAT0=A    W        ;STORE ABSVAL IN ANSWER ARRAY
  849.          D0=D0+    16       ;POINT AT NEXT ANSWER
  850. * CHECK FOR DONE
  851.          C=R4.F    A        ;GET COUNTER
  852.          SETHEX
  853.          C=C-1     A
  854.          R4=C.F    A        ;DEC COUNTER
  855.          GONC    ABSLOOP    ;LOOP UNTIL BORROW
  856.  
  857. ***NOW ANSWER CONTAINS ABS VALUES OF ORIG
  858.          GOVLNG    =GPOverWrR0Lp      ;GETPTR, REPLACE TOS WITH ANS, LOOP
  859.  
  860. ************************ SQDAT1 ***********************************
  861. * GET REAL FROM DAT1, UNPACK AND SQUARE IT, INC D1
  862.  
  863. SQDAT1   A=DAT1    W        ;A=REAL
  864.          D1=D1+    16       ;D1->NEXT
  865.          GOSBVL    =SPLITA  ;AB=REAL
  866.          C=B       W
  867.          D=C       W
  868.          C=A       W        ;COPY INTO CD
  869.          GOSBVL    =MULTF   ;AB=SQR(REAL)
  870.          RTN
  871.  
  872. ENDCODE
  873. ;
  874. ;
  875.